function beta_p_mat = sim_init_betap(data,fp,col_long,betap_draws)

% Generate random draws of beta_p for every simulated household (N households, R draws per household)
% Draws are correlated with parental education level (taken as exogenous)

beta_p_mat = NaN(fp.N,fp.R);

for i = 1:fp.N
    % household i's data
    data_here = data(1+(fp.M+1)*(i-1):(fp.M+1)*i, : );

    if fp.compstat_SES_homprefs == 1
        % In this comparative exercise, we shut down heterogeneity in parental discount factor distribution by SES
        % Note: we still allow for unconditional heterogeneity (iid random draws from the unconditional distribution).
        dad_educ = fp.s2_median;
    elseif fp.compstat_SES_homprefs == 0
        dad_educ = data_here(1,col_long.dad_educ);
    end

    % First determine parental education level, for beta_p
    % This assumes fp.N_educ_betap = 3 categories
    if dad_educ <= 12
        jj = 1; % high school or less
    elseif dad_educ >= 13 && dad_educ <= 16
        jj = 2; % some college or college
    elseif dad_educ >= 17
        jj = 3; % graduate
    end

    for r = 1:fp.R
        % Randomly draw discount factor for every (q,r), keeping education level fixed across all r draws
        beta_p_mat(i,r) = fp.betap_vec(find(betap_draws(i,r)<=fp.betap_cdf(:,jj),1));
    end
end

end % function